#include <cmath>
#include <sys/types.h>
#include <sys/stat.h>
#include <fstream>
#include <limits>
#include <cstdio>
#include <cassert>
#include <unistd.h>
#include <cerrno>
#include "bdbibd.h"
#include "xstdlib.h"
#include "xzlib.h"
#include "where.h"
#include "xstring.h"

using namespace std;
using namespace bdb;



#define FAMILY_DB_NAME "id2family.db"
#define INDIVIDUAL_DB_NAME "id2individual.db"
#define CHROMOSOME_DB_NAME "id2chromosome.db"
#define PAIROFINDI_TO_PAIRID_NAME "indi2pairid.db"
#define PAIRID_TO_PAIROFINDI_NAME "pairid2indi.db"
#define PAIR_ID_TO_THEORITICAL "pairId2theoritical.db"
#define LINKAGE_DB_NAME "linkage.db"
#define RS_TO_POSITION "rs2position.db"
#define POSITION_TO_MARKER "position2marker.db"
#define METADATADB "metadata.db"

#define COL_PROBSET_ID "Probe Set ID"
#define COL_DBSNP_RS_ID "dbSNP RS ID"
#define COL_CHROMOSOME "Chromosome"
#define COL_POSITION "Physical Position"
#define COL_ALLELE_A "Allele A"
#define COL_ALLELE_B "Allele B"
#define COL_GENETIC_MAP "Genetic Map"

#define LOG(a) do { if(debug) {std::cerr << __FILE__ <<":"<<__FUNCTION__ <<"["<< __LINE__ << "]" << a << endl;}} while(false)




std::ostream& operator<<(std::ostream& out,const PairOfIndividual& o)
    {
    out << "pair("<< o.first<<","<< o.second <<")";
    return out;
    }

std::ostream& operator<<(std::ostream& out,const Genotype& o)
    {
    for(uint32_t i=0;i< COUNT_PVALUE;++i)
	{
	if(i>0) out << "|";
	out << o.pvalue[i];
	}
    return out;
    }

std::ostream& operator<<(std::ostream& out,const Position& o)
    {
    out << (int)o.chrom << ":" << o.pos;
    return out;
    }

std::ostream& operator<<(std::ostream& out,const Marker& o)
    {
    out << o.rsId << ";" << o.position << ";";
    for(size_t i=0;i< o.geneticMaps.size();++i)
	{
	if(i>0) out << "|";
	out << o.geneticMaps.at(i);
	}
    return out;
    }

static void split(vector<char*>& v,char* str,const char* delim)
	{
	v.clear();
	char* p2=NULL;
	char *scan=strtok_r(str,delim, &p2);
	while(scan!=NULL)
		{
		v.push_back(scan);
		scan=strtok_r(NULL,delim, &p2);
		}
	}

/***********************************************************************************/
/***********************************************************************************/
/***********************************************************************************/
CoupleOfIndividual::CoupleOfIndividual():
		first((IndividualPtr)NULL),
		second((IndividualPtr)NULL)
	{
	}

CoupleOfIndividual::CoupleOfIndividual(IndividualPtr i1,IndividualPtr i2)
	{
	if(i1==NULL || i2==NULL) THROW("Null ptr exception");
	if(i1->id< i2->id)
		{
		first=i1;
		second=i2;
		}
	else
		{
		first=i2;
		second=i1;
		}
	}

bool CoupleOfIndividual::contains(IndividualPtr id) const
	{
	return first==id || second==id;
	}

IndividualPtr CoupleOfIndividual::complementary(IndividualPtr id) const
	{
	if(first==id) return second;
	if(second==id) return first;
	return NULL;
	}

int32_t CoupleOfIndividual::concordant() const
	{
	if(first->status==2 && second->status==2)
		{
		return 2;
		}
	else if((first->status==1 && second->status==2) ||
			(first->status==2 && second->status==1))
		{
		return 1;
		}
	else
		{
		return 0;
		}
	}


/***********************************************************************************/
/***********************************************************************************/
/***********************************************************************************/
Genotype::Genotype()
    {
    for(int i=0;i< COUNT_PVALUE;++i)
	{
	pvalue[i]=-1;
	}
    }

bool Genotype::hasData() const
	{
	for(int i=0;i< COUNT_PVALUE;++i)
		{
		if(pvalue[i]>=0) return true;
		}
	return false;
	}

pvalue_t Genotype::ibd0() const
	{
	return 	ibd(IBD_0);
	}

pvalue_t Genotype::ibd(WhichIBD w) const
	{
	switch(w)
	    {
	    case IBD_0:
	    case IBD_1:
	    case IBD_2: return pvalue[w];
	    default:THROW("ERROR");
	    }
	THROW("ERROR");
	}


/***********************************************************************************/
/***********************************************************************************/
/***********************************************************************************/

PairOfIndividual::PairOfIndividual():first((individual_id_t)-1),
	second((individual_id_t)-1)
    {
    }

PairOfIndividual::PairOfIndividual(individual_id_t i1,individual_id_t i2)
	{
	if(i1< i2)
		{
		first=i1;
		second=i2;
		}
	else
		{
		first=i2;
		second=i1;
		}
	}

bool PairOfIndividual::contains(individual_id_t id) const
	{
	return first==id || second==id;
	}

individual_id_t PairOfIndividual::complementary(individual_id_t id) const
	{
	if(first==id) return second;
	if(second==id) return first;
	return (individual_id_t)-1;
	}
/***********************************************************************************/
/***********************************************************************************/
/***********************************************************************************/


LinkageRow::LinkageRow(uint32_t size):len(size),array(0)
    {
    array=new Genotype[len];
    Genotype g;
    for(int i=0;i< COUNT_PVALUE;++i)
	{
	g[i]=-1;
	}
    for(uint32_t i=0;i< len;++i)
	{
	array[i]=g;
	}
    }

LinkageRow::LinkageRow(const bdb::Entry* entry):len(0)
    {
    if(entry->size()%sizeof(Genotype)!=0) THROW("Error");
    len=entry->size()/sizeof(Genotype);
    array=new Genotype[len];
    memcpy((void*)array,entry->data(),entry->size());
    }


LinkageRow::~LinkageRow()
    {
    if(array!=0) delete[] array;
    }

const Genotype* LinkageRow::getGenotype(pair_id_t pairid) const
    {
    if(pairid>=size()) THROW("out of bound");
    return &array[pairid];
    }

pvalue_t LinkageRow::ibd(pair_id_t pairid,WhichIBD whichibd) const
    {
    return getGenotype(pairid)->ibd(whichibd);
    }

pvalue_t LinkageRow::ibd0(pair_id_t pairid) const
    {
    return ibd(pairid,IBD_0);
    }


void LinkageRow::setGenotype(pair_id_t pairid,const Genotype* genotype)
    {
    if(size()< pairid)
	    {
	    Genotype* array2=new Genotype[pair_id+1];
	    std::copy(&array[0],&array[len],&array2[0]);
	    uint32_t oldsize=size();
	    delete [] array;
	    this->array=array2;
	    this->len=pairid+1;
	    Genotype nil;
	    for(int i=oldsize;i< pair_id;++i)
		    {
		    setGenotype(i,nil);
		    }
	    }
    array[pairid]=*(genotype);
    }


/***********************************************************************************/
/***********************************************************************************/
/***********************************************************************************/

IBDEnv::IBDEnv():txn(NULL),
	familyDb(NULL),
	individualsDb(NULL),
	chromName2TidDb(NULL),
	pairIdToPairOfIndividualsDb(NULL),
	pairOfIndividualsToPairId(NULL),
	linkageDB(NULL),
	rsToPositionDb(NULL),
	pos2markerDb(NULL),
	pairId2theoriticalDb(NULL),
	metadataDB(NULL),
	count_pairs(0U),
	debug(false)
	{

	}
#define CLOSE_BDB(a) LOG(#a); if(a!=NULL) { LOG("closing" << #a); a->close(); LOG("done " << #a); a->close(); delete a; a=NULL;}

void IBDEnv::close()
	{
	for(map<family_id_t,Family*>::iterator r=id2family.begin();
		r!=id2family.end();
		++r)
		{
		delete r->second;
		}
	id2family.clear();
	name2family.clear();
	for(map<individual_id_t,Individual*>::iterator r=id2individual.begin();
		r!=id2individual.end();
		++r)
		{
		delete r->second;
		}

	id2individual.clear();
	id2chromosome.clear();
	chrom2tid.clear();
	count_pairs=0;

	CLOSE_BDB(rsToPositionDb);
	CLOSE_BDB(pos2markerDb);
	CLOSE_BDB(familyDb);
	CLOSE_BDB(individualsDb);
	CLOSE_BDB(chromName2TidDb);
	CLOSE_BDB(pairIdToPairOfIndividualsDb);
	CLOSE_BDB(pairId2theoriticalDb);
	CLOSE_BDB(pairOfIndividualsToPairId);
	CLOSE_BDB(linkageDB);
	CLOSE_BDB(metadataDB);
	env.close();
	}

#undef CLOSE_BDB

IBDEnv::~IBDEnv()
	{
	close();
	}

Family* IBDEnv::entryToFamily(const Entry* entry)
	{
	EntryReader r(entry);
	family_id_t id;
	r.read((void*)&id,sizeof(family_id_t));
	std::string name=r.readString8();
	return new Family(id,name);
	}

void IBDEnv::familyToEntry(const FamilyPtr family,Entry* entry)
	{
	EntryWriter w(entry);
	w.write((void*)&family->id,sizeof(family_id_t));
	w.writeString8(family->name);
	}

Individual* IBDEnv::entryToIndividual(const Entry* entry)
	{
	EntryReader r(entry);
	family_id_t family_id;
	individual_id_t id;
	r.read((void*)&id,sizeof(individual_id_t));
	r.read((void*)&family_id,sizeof(family_id_t));
	map<family_id_t,Family*>::iterator iter=id2family.find(family_id);
	if(iter==id2family.end()) THROW("Cannot find family id " << family_id);
	Family* fam=iter->second;
	string name=r.readString8();
	string fatherId=r.readString8();
	string motherId=r.readString8();
	Gender gender;
	r.read((void*)&gender,sizeof(Gender));
	phenotype_t phenotype;
	r.read((void*)&phenotype,sizeof(phenotype_t));
	return new Individual(id,fam,name,fatherId,motherId,gender,phenotype);
	}

void IBDEnv::individualToEntry(const IndividualPtr indi,Entry* entry)
	{
	EntryWriter w(entry);
	w.write((void*)&(indi->id),sizeof(individual_id_t));
	w.write((void*)&(indi->family->id),sizeof(family_id_t));
	w.writeString8(indi->name);
	w.writeString8(indi->father);
	w.writeString8(indi->mother);
	w.write(&(indi->gender),sizeof(Gender));
	w.write(&(indi->status),sizeof(phenotype_t));
	}


Marker* IBDEnv::entryToMarker(const Entry* entry)
	{
	EntryReader r(entry);
	Marker* marker=new Marker;
	marker->rsId=r.readString8();
	r.read((void*)&(marker->position),sizeof(Position));
	uint8_t n=r.readUInt8();
	for(int i=0;i< n;++i)
		{
		genetic_map_t gm;
		r.read((void*)&gm,sizeof(genetic_map_t));
		marker->geneticMaps.push_back(gm);
		}
	return marker;
	}

void IBDEnv::markerToEntry(const Marker* marker,Entry* entry)
	{
	EntryWriter w(entry);
	w.writeString8(marker->rsId);
	w.write(&(marker->position),sizeof(Position));
	w.writeUInt8((uint8_t)marker->geneticMaps.size());
	for(size_t i=0;i< marker->geneticMaps.size();++i)
		{
		genetic_map_t gm=marker->geneticMaps.at(i);
		w.write((void*)&gm,sizeof(genetic_map_t));
		}
	}

Family* IBDEnv::findFamilyByName(const char* familyName)
		{
		if(familyName==NULL || familyName[0]==0) return NULL;
		map<std::string,Family*>::iterator r=name2family.find(familyName);
		if(r==name2family.end()) return NULL;
		return r->second;
		}

int IBDEnv::findTidByName(const char* chromName)
	{
	if(chromName==NULL || chromName[0]==0) return -1;
	map<std::string,tid_t>::iterator r=chrom2tid.find(chromName);
	if(r==chrom2tid.end()) return -1;
	return r->second;
	}

int IBDEnv::comparePositions(DB *db, const DBT *dbt1, const DBT *dbt2)
	{
	Position p1;
	Position p2;
	assert(dbt1->size>=sizeof(Position));
	assert(dbt2->size>=sizeof(Position));
	memcpy(&p1,(void*)dbt1->data,sizeof(Position));
	memcpy(&p2,(void*)dbt2->data,sizeof(Position));
	int i= p1.chrom -p2.chrom;
	if(i!=0) return i;
	return p1.pos-p2.pos;
	}

int IBDEnv::comparePairs(DB *db, const DBT *dbt1, const DBT *dbt2)
	{
	PairOfIndividual p1;
	PairOfIndividual p2;
	assert(dbt1->size>=sizeof(PairOfIndividual));
	assert(dbt2->size>=sizeof(PairOfIndividual));
	memcpy(&p1,(void*)dbt1->data,sizeof(PairOfIndividual));
	memcpy(&p2,(void*)dbt2->data,sizeof(PairOfIndividual));
	int i= p1.first -p2.first;
	if(i!=0) return i;
	return p1.second-p2.second;
	}


IndividualPtr IBDEnv::findIndividualById(individual_id_t id)
    {
    map<individual_id_t,Individual*>::iterator ri=this->id2individual.find(id);
    return(ri==this->id2individual.end()?NULL:ri->second);
    }

CoupleOfIndividual IBDEnv::findIndividualsByPair(PairOfIndividual* pair)
    {
    return CoupleOfIndividual(
	    findIndividualById(pair->first),
	    findIndividualById(pair->second)
	    );
    }

pair_id_t IBDEnv::findPairIdByPair(PairOfIndividual* pair)
    {
	Entry k((void*)pair,sizeof(PairOfIndividual));
    auto_ptr<Entry> p=pairOfIndividualsToPairId->get(txn,&k);
    if(p.get()==NULL) return (pair_id_t)-1;
    pair_id_t id;
    assert(p->size()==sizeof(pair_id_t));
    memcpy((void*)&id,(void*)p->data(),sizeof(pair_id_t));
    return id;
    }


auto_ptr<Marker> IBDEnv::findMarkerByPosition(const Position* position)
    {
    Entry posEntry((void*)position,sizeof(Position));
    auto_ptr<Entry> markerEntry= this->pos2markerDb->get(this->txn,&posEntry,0);
    if(markerEntry.get()==NULL)
		{
   		return auto_ptr<Marker>();
		}
    return auto_ptr<Marker>(this->entryToMarker(markerEntry.get()));
    }


auto_ptr<LinkageRow> IBDEnv::getLinkageRowByMarker(const char* rs_name)
    {
//TODO
    }

std::auto_ptr<Genotype> IBDEnv::findTheoriticalIBDByPairId(pair_id_t id)
	{
	Entry k((void*)&id,sizeof(pair_id_t));
	auto_ptr<Entry> e= this->pairId2theoriticalDb->get(this->txn,&k,0);
	if(e.get()==NULL)
		{
		return auto_ptr<Genotype>();
		}
	Genotype* g=new Genotype;
	assert(e->size()>= sizeof(Genotype));
	memcpy((void*)g,e->data(),sizeof(Genotype));
	return auto_ptr<Genotype>(g);
	}


void IBDEnv::open(char* dir,bool readOnly)
	{
	this->readOnly=readOnly;
	if(!readOnly)
			{
			//only effects newly created environments
			env.setCacheSize(0,64 * 1024 * 1024);//64Mo
			}
	env.open(dir,(readOnly?0:DB_CREATE)|DB_INIT_MPOOL,0);

	familyDb=new Database(env.env());
	familyDb->setBTreeCompare(db_sorter<family_id_t>::comparator());
	individualsDb=new Database(env.env());
	individualsDb->setBTreeCompare(db_sorter<individual_id_t>::comparator());
	chromName2TidDb=new Database(env.env());
	pairOfIndividualsToPairId=new Database(env.env());
	pairOfIndividualsToPairId->setBTreeCompare(IBDEnv::comparePairs);
	pairId2theoriticalDb=new Database(env.env());
	pairId2theoriticalDb->setBTreeCompare(db_sorter<pair_id_t>::comparator());
	pairIdToPairOfIndividualsDb=new Database(env.env());
	pairIdToPairOfIndividualsDb->setBTreeCompare(db_sorter<pair_id_t>::comparator());
	linkageDB=new Database(env.env());
	linkageDB->setBTreeCompare(IBDEnv::comparePositions);

	rsToPositionDb=new Database(env.env());
	pos2markerDb=new Database(env.env());


	pos2markerDb->setBTreeCompare(IBDEnv::comparePositions);
	if(!readOnly)
	    {
#define DEFAULT_OPEN_DB_WRITE(var,name)	if(var==NULL) THROW("Illegal State"); var->open(txn,name,name, DB_BTREE, DB_CREATE|DB_TRUNCATE);
	    DEFAULT_OPEN_DB_WRITE(familyDb,FAMILY_DB_NAME)
	    DEFAULT_OPEN_DB_WRITE(individualsDb,INDIVIDUAL_DB_NAME)
	    DEFAULT_OPEN_DB_WRITE(pos2markerDb,POSITION_TO_MARKER)
	    DEFAULT_OPEN_DB_WRITE(rsToPositionDb,RS_TO_POSITION)
	    DEFAULT_OPEN_DB_WRITE(chromName2TidDb,CHROMOSOME_DB_NAME)
	    DEFAULT_OPEN_DB_WRITE(pairIdToPairOfIndividualsDb,PAIRID_TO_PAIROFINDI_NAME)
	    DEFAULT_OPEN_DB_WRITE(pairId2theoriticalDb,PAIR_ID_TO_THEORITICAL)
	    DEFAULT_OPEN_DB_WRITE(pairOfIndividualsToPairId,PAIROFINDI_TO_PAIRID_NAME)
	    DEFAULT_OPEN_DB_WRITE(linkageDB,LINKAGE_DB_NAME)
	    DEFAULT_OPEN_DB_WRITE(metadataDB,METADATADB)
#undef DEFAULT_OPEN_DB_WRITE
	    }

	else /* read only */
		{
#define DEFAULT_OPEN_DB_RDONLY(var,name) if(var==NULL) THROW("Illegal State"); var->open(txn,name,name, DB_UNKNOWN, DB_RDONLY)
		DEFAULT_OPEN_DB_RDONLY(familyDb,FAMILY_DB_NAME);
		/* family */
		std::auto_ptr<Cursor> c= familyDb->cursor(txn);
		while(c->next())
			{
			Family* fam= this->entryToFamily(c->dataEntry().get());
			id2family.insert(make_pair(fam->id,fam));
			name2family.insert(make_pair(fam->name,fam));
			}
		c.reset();
		familyDb->close();


		/* individuals */
		DEFAULT_OPEN_DB_RDONLY(individualsDb,INDIVIDUAL_DB_NAME);
		c= individualsDb->cursor(txn);
		while(c->next())
			{
			Individual* indi= entryToIndividual(c->dataEntry().get());
			id2individual.insert(make_pair(indi->id,indi));
			indi->family->put(indi);
			}
		c.reset();
		individualsDb->close();


		/* chromosomes */
		DEFAULT_OPEN_DB_RDONLY(chromName2TidDb,CHROMOSOME_DB_NAME);
		c= chromName2TidDb->cursor(txn);
		while(c->next())
			{
			tid_t cid;
			string name=c->keyToString();
			memcpy((void*)&cid,c->data()->data,sizeof(tid_t));
			id2chromosome.insert(make_pair(cid,name));
			chrom2tid.insert(make_pair(name,cid));
			}
		c.reset();
		chromName2TidDb->close();

		/* pairs */
		DEFAULT_OPEN_DB_RDONLY(pairId2theoriticalDb,PAIR_ID_TO_THEORITICAL);
		DEFAULT_OPEN_DB_RDONLY(pairIdToPairOfIndividualsDb,PAIRID_TO_PAIROFINDI_NAME);
		this->count_pairs=(uint32_t)pairIdToPairOfIndividualsDb->count(txn,0);
		/*open other DB */
		DEFAULT_OPEN_DB_RDONLY(rsToPositionDb,RS_TO_POSITION);
		DEFAULT_OPEN_DB_RDONLY(pos2markerDb,POSITION_TO_MARKER);
		DEFAULT_OPEN_DB_RDONLY(pairOfIndividualsToPairId,PAIROFINDI_TO_PAIRID_NAME);
		DEFAULT_OPEN_DB_RDONLY(linkageDB,LINKAGE_DB_NAME);
		DEFAULT_OPEN_DB_RDONLY(metadataDB,METADATADB);
#undef DEFAULT_OPEN_DB_RDONLY
		}
	
	}


void IBDEnv::putMetadata(std::string key,std::string value)
    {
    Entry keyE(key);
    Entry dataE(value);
    metadataDB->put(txn,&keyE,&dataE,0);
    }


void IBDEnv::mappingSplit(const std::string& s,std::vector<std::string>& tokens)
	{
	tokens.clear();
	char delim='\0';
	std::string::size_type i=0;
	while(i< s.size())
		{
		if(!(s[i]=='\'' || s[i]=='\"')) THROW("expected a quote or apos in "<<s);
		delim=s[i];
		i++;
		if(i==s.size())
			{
			THROW("unfinished quoted string " << s);
			}
		std::string::size_type beg=i;
		std::string::size_type end=i;
		for(;;)
		    {
		    if(i==s.size())
			{
				THROW("unfinished quoted string " << s);
			}
		     if(s[i]==delim)
			{
				end=i;
				i++;
				break;
			}
			i++;
		    }
		std::string tok(s.substr(beg,end-beg));
		//LOG("tokens["<< tokens.size()<< "]=" << tok );
		tokens.push_back(tok);
		if(i==s.size()) break;
		if(s[i]!=',') THROW("comma missing in "<< s);
		i++;
		if(!(s[i]=='\'' || s[i]=='\"')) THROW("expected a quote at the end of "<< s);
		}
	}


void IBDEnv::_prefillTids()
    {
    for(tid_t i=1;i< 23;++i)
	{
	char tmp[10];
	sprintf(tmp,"%d",i);
	chrom2tid.insert(make_pair(tmp,i));
	sprintf(tmp,"chr%d",i);
	this->id2chromosome.insert(make_pair(i,tmp));
	}
    chrom2tid.insert(make_pair("X",23));
    chrom2tid.insert(make_pair("chrX",23));
    this->id2chromosome.insert(make_pair(23,"chrX"));
    chrom2tid.insert(make_pair("Y",24));
    chrom2tid.insert(make_pair("chrY",24));
    this->id2chromosome.insert(make_pair(24,"chrX"));
    }

#define SCRAP_HEADER(themacro,var)  if(header.at(i).compare(themacro)==0) \
	{\
	if(var!=-1) THROW("Duplicate "<<themacro<<" in "<<line);\
	var=i;\
	}
#define CHECK_HEADER(var,themacro)  if(var==-1) THROW("Cannot get " << themacro << " in " << line)
void IBDEnv::readMapping(const char* file)
	{
	size_t n_insert=0;
	int colIndexRsId=-1;
	int colChrom=-1;
	int colPosition=-1;
	int colProbSetId=-1;
	int colAlleleA=-1;
	int colAlleleB=-1;
	int colGeneticMap=-1;
	LOG( "Reading mapping file " << file );

	_prefillTids();

	vector<string> header;
	GZReader in(file);
	char* row=NULL;


	while((row=in.readLine())!=NULL)
		{
		if(row[0]==0 || row[0]=='#') continue;
		string line(row);
		mappingSplit(line,header);

		for(size_t i=0;i< header.size();++i)
			{
			SCRAP_HEADER(COL_CHROMOSOME,colChrom)
			else SCRAP_HEADER(COL_POSITION,colPosition)
			else SCRAP_HEADER(COL_DBSNP_RS_ID,colIndexRsId)
			else SCRAP_HEADER(COL_PROBSET_ID,colProbSetId)
			else SCRAP_HEADER(COL_ALLELE_A,colAlleleA)
			else SCRAP_HEADER(COL_ALLELE_B,colAlleleB)
			else SCRAP_HEADER(COL_GENETIC_MAP,colGeneticMap)
			}
		CHECK_HEADER(colChrom,COL_CHROMOSOME);
		CHECK_HEADER(colPosition,COL_POSITION);
		CHECK_HEADER(colIndexRsId,COL_DBSNP_RS_ID);
		CHECK_HEADER(colProbSetId,COL_PROBSET_ID);
		CHECK_HEADER(colAlleleA,COL_ALLELE_A);
		CHECK_HEADER(colAlleleB,COL_ALLELE_B);
		CHECK_HEADER(colGeneticMap,COL_GENETIC_MAP);
		break;
		}
	if(header.empty())
		{
		THROW("header not found in "<< file);
		}
	std::vector<std::string> tokens;
	while((row=in.readLine())!=NULL)
		{
		string line(row);
		mappingSplit(line,tokens);
		if(tokens.size()!=header.size())
			{
			THROW("expected "<< header.size()
				<< " columns but got "
				<< tokens.size() <<" in "+line);
			}
		if(tokens[colChrom].compare("---")==0 ||
		   tokens[colPosition].compare("---")==0)
			  {
			  //ignore
			  continue;
			  }

		Marker* marker=new Marker;

		/* chromosome */
		map<string,tid_t>::iterator rc=chrom2tid.find(tokens.at(colChrom));
		if(rc==chrom2tid.end())
			{
			marker->position.chrom=(tid_t)(1+chrom2tid.size());
			if(id2chromosome.find(marker->position.chrom)!=id2chromosome.end()) THROW("Runtime error");
			this->chrom2tid.insert(make_pair(tokens.at(colChrom),marker->position.chrom));
			this->id2chromosome.insert(make_pair(marker->position.chrom,tokens.at(colChrom)));
			}
		else
			{
			marker->position.chrom=rc->second;
			}


		marker->rsId=tokens.at(colIndexRsId);



		marker->position.pos=(position_t)std::atol(tokens.at(colPosition).c_str());
		//genetic map
		string::size_type tripleslash=0;
		while(tripleslash<tokens[colGeneticMap].size())
			{
			/* example: token to parser:
			26.2497537844838 // D8S1827 // D8S549 // --- // --- /// 31.1074557729995 // D8S1827 // D8S549 // AFM107YA1 // AFM303ZC1 /// 28.7621103124195 // --- // --- // 259984 // 537009
			*/

			std::string::size_type doubeslash=tokens[colGeneticMap].find(" // ",tripleslash);
			if(doubeslash==std::string::npos)
				{
				break;
				}
			marker->geneticMaps.push_back((genetic_map_t)atof(tokens[colGeneticMap].substr(tripleslash,doubeslash-tripleslash).c_str()));
			tripleslash=doubeslash;


			tripleslash = tokens[colGeneticMap].find(" /// ",tripleslash);
			if(tripleslash==std::string::npos)
				{
				break;
				}
			tripleslash+=5;
			}

		//insert pos2markerDb
			{
			Entry keyE((void*)&(marker->position),sizeof(Position));
			if(pos2markerDb->contains(txn,&keyE,0))
				{
				LOG( "duplicate pos/rs "<< line );
				}
			Entry dataE;
			markerToEntry(marker,&dataE);
			pos2markerDb->put(txn,&keyE,&dataE,0);
			//if(!pos2markerDb->contains(txn,&keyE,0)) THROW("BOUM");
			}
		//insert rs2pos
			{
			Entry keyE(marker->rsId.c_str());
			if(rsToPositionDb->contains(txn,&keyE,0))
				{
				LOG( "duplicate rs/position"<< line );
				}
			Entry dataE((void*)&(marker->position),sizeof(Position));
			rsToPositionDb->put(txn,&keyE,&dataE,0);
			//if(!rsToPositionDb->contains(txn,&keyE,0)) THROW("BOUM");
			}
		n_insert++;
		if(n_insert%10000==0)
			{
			LOG("Marker inserted:"<< n_insert << " current:" << marker->rsId );
			}
		delete marker;
		}
	in.close();
	rsToPositionDb->sync(0);

#ifdef PARANOID
	auto_ptr<Cursor> c= rsToPositionDb->cursor();
	while(c->next())
		{
		LOG(c->keyToString());
		Position p;
		memcpy(&p,c->data()->data,sizeof(Position));
		LOG(p);
		}
	c.reset();

	c= pos2markerDb->cursor();
	while(c->next())
			{
			Position p;
			memcpy(&p,c->key()->data,sizeof(Position));
			LOG(p);
			Marker* m=entryToMarker(&(*c->dataEntry()));
			LOG(*m);
			delete m;
			}
	c.reset();
#endif

	//save chrom2tid
	for(map<string,tid_t>::iterator rc=this->chrom2tid.begin();
		rc!=chrom2tid.end();
		++rc)
		{
		Entry keyE((void*)rc->first.data(),rc->first.size());
		Entry dataE((void*)&(rc->second),sizeof(tid_t));

		this->chromName2TidDb->put(txn,&keyE,&dataE);
		}
	chromName2TidDb->sync(0);
	}



void IBDEnv::readBed(const char* file)
	{
	size_t n_insert=0;

	const int colChrom=0;
	const int colStart=1;
	const int colEnd=2;
	const int colIndexRsId=3;
	LOG( "Reading BED file " << file );

	_prefillTids();

	vector<string> header;
	GZReader in(file);
	char* row=NULL;


	std::vector<std::string> tokens;
	while((row=in.readLine())!=NULL)
		{
		string line(row);
		if(line.empty() || line[0]=='#') continue;

		split(tokens,line,"\t");
		if(tokens.size()<4)
		    {
		    THROW("BED: illegal number of tokens in "<< line);
		    }

		if(tokens[0].empty())
		    {
		    THROW("BED: bad chromosome "<< line);
		    }
		if(tokens[1].empty() || tokens[2].empty())
		    {
		    THROW("BED: bad start/end "<< line);
		    }

		this->

		Marker* marker=new Marker;

		/* chromosome */
		map<string,tid_t>::iterator rc=chrom2tid.find(tokens.at(colChrom));
		if(rc==chrom2tid.end())
			{
			marker->position.chrom=(tid_t)(1+chrom2tid.size());
			if(id2chromosome.find(marker->position.chrom)!=id2chromosome.end()) THROW("Runtime error");
			this->chrom2tid.insert(make_pair(tokens.at(colChrom),marker->position.chrom));
			this->id2chromosome.insert(make_pair(marker->position.chrom,tokens.at(colChrom)));
			}
		else
			{
			marker->position.chrom=rc->second;
			}


		marker->rsId=tokens.at(colIndexRsId);
		marker->position.pos=(position_t)std::atol(tokens.at(colStart).c_str());

		//insert pos2markerDb
			{
			Entry keyE((void*)&(marker->position),sizeof(Position));
			if(pos2markerDb->contains(txn,&keyE,0))
				{
				LOG( "duplicate pos/rs "<< line );
				}
			Entry dataE;
			markerToEntry(marker,&dataE);
			pos2markerDb->put(txn,&keyE,&dataE,0);
			//if(!pos2markerDb->contains(txn,&keyE,0)) THROW("BOUM");
			}
		//insert rs2pos
			{
			Entry keyE(marker->rsId.c_str());
			if(rsToPositionDb->contains(txn,&keyE,0))
				{
				LOG( "duplicate rs/position"<< line );
				}
			Entry dataE((void*)&(marker->position),sizeof(Position));
			rsToPositionDb->put(txn,&keyE,&dataE,0);
			//if(!rsToPositionDb->contains(txn,&keyE,0)) THROW("BOUM");
			}
		n_insert++;
		if(n_insert%10000==0)
			{
			LOG("Marker inserted:"<< n_insert << " current:" << marker->rsId );
			}
		delete marker;
		}
	in.close();
	rsToPositionDb->sync(0);

#ifdef PARANOID
	auto_ptr<Cursor> c= rsToPositionDb->cursor();
	while(c->next())
		{
		LOG(c->keyToString());
		Position p;
		memcpy(&p,c->data()->data,sizeof(Position));
		LOG(p);
		}
	c.reset();

	c= pos2markerDb->cursor();
	while(c->next())
			{
			Position p;
			memcpy(&p,c->key()->data,sizeof(Position));
			LOG(p);
			Marker* m=entryToMarker(&(*c->dataEntry()));
			LOG(*m);
			delete m;
			}
	c.reset();
#endif

	//save chrom2tid
	for(map<string,tid_t>::iterator rc=this->chrom2tid.begin();
		rc!=chrom2tid.end();
		++rc)
		{
		Entry keyE((void*)rc->first.data(),rc->first.size());
		Entry dataE((void*)&(rc->second),sizeof(tid_t));

		this->chromName2TidDb->put(txn,&keyE,&dataE);
		}
	chromName2TidDb->sync(0);
	}



void IBDEnv::readPedigree(const char* file)
	{
	vector<char*> tokens;
	GZReader in(file);
	char* line=NULL;

	while((line=in.readLine())!=NULL)
		{
		if(line[0]==0 || line[0]=='#') continue;

		split(tokens,line," \t");
		if(tokens.size()!=6) THROW("illegal number of tokens in "<< line);
		Family* family=NULL;
		map<string,Family*>::iterator r=name2family.find(tokens[0]);
		if(r==name2family.end())
			{
			family=new Family(name2family.size()+1,tokens[0]);
			name2family.insert(make_pair(family->name,family));
			id2family.insert(make_pair(family->id,family));
			}
		else
			{
			family=r->second;
			}
		string name=tokens[1];
		string fatherId=tokens[2];
		string motherId=tokens[3];
		int gender=atoi(tokens[4]);
		switch((Gender)gender)
			{
			case gender_male:
			case gender_female:
			case gender_unknown:break;
			default:THROW("unknow gender "<< tokens[4]);
			}
		phenotype_t status=atoi(tokens[5]);
		Individual* indi=new Individual(
				id2individual.size()+1,
				family,
				name,
				fatherId,motherId,
				(Gender)gender,status);
		id2individual.insert(make_pair(indi->id,indi));
		family->put(indi);

		}
	in.close();
	/* Family */

	for(map<family_id_t,Family*>::iterator r=id2family.begin();
		r!=id2family.end();
		++r)
		{
		Entry key;
		Entry data;
		key.write((void*)&(r->second->id),sizeof(family_id_t));
		familyToEntry(r->second,&data);
		familyDb->put(txn,&key,&data);
		}
	familyDb->sync();
	
	

	/* Individual */
	for(map<individual_id_t,Individual*>::iterator r=id2individual.begin();
		r!=id2individual.end();
		++r)
		{
		Entry key;
		Entry data;
		key.write((void*)&(r->second->id),sizeof(individual_id_t));
		individualToEntry(r->second,&data);
		individualsDb->put(txn,&key,&data);
		}
	individualsDb->sync();
	}


void IBDEnv::readReskin(const char* file)
	{

	GZReader in(file);
	char* line;
	while((line=in.readLine())!=NULL)
			{
			if(line[0]==0 || line[0]=='#') continue;
			vector<string> tokens;
			XCString::split(line,"\t",tokens);
			if(tokens.size()!=12) THROW("Expected 12 columns in "<< line);
			Family* fam=findFamilyByName(tokens[0].c_str());
			if(fam==NULL) THROW("Cannot find family "<< tokens[0]);
			Individual* indi1=fam->findIndividualByName(tokens[1]);
			if(indi1==NULL) THROW("Cannot find individual1 "<< tokens[1]);
			Individual* indi2=fam->findIndividualByName(tokens[2]);
			if(indi2==NULL) THROW("Cannot find individual2 "<< tokens[2]);
			PairOfIndividual couple(indi1->id,indi2->id);
			Genotype theoritical;
			memset((void*)&theoritical,0,sizeof(Genotype));
			for(int i=0;i< COUNT_PVALUE;++i)
				{
				char* p2;
				theoritical.pvalue[i]=strtod(tokens[11-i].c_str(),&p2);//inverse (fichier is ibd 2 ; 1; 0 )
				if(*p2!=0) THROW("Bad number "<< tokens[9+i]);
				}

			pair_id_t pair_id=findPairIdByPair(&couple);
			if(pair_id!=(pair_id_t)-1)//only if we saw this pair in 'getData'
				{
				Entry dataE((void*)&pair_id,sizeof(pair_id_t));
				Entry theoE((void*)&theoritical,sizeof(Genotype));
				//add theo data
				pairId2theoriticalDb->put(txn,&dataE,&theoE,0);
				}
			}
	in.close();
	pairId2theoriticalDb->sync();
	}


#ifdef BAD_VERSION
void IBDEnv::readData(std::set<std::string>& datafiles)
	{
	int max_warning_count=10;
	/** row of count_pair*Genotype */
	size_t size_genotype_row=count_pairs*sizeof(Genotype);;
	Genotype* genotypeRow=(Genotype*)safeCalloc(count_pairs,sizeof(Genotype));

	string prev_fam;
	Family* prev_fam_ptr=NULL;
	string prev_indi1;
	Individual* prev_indi1_ptr=NULL;
	string prev_indi2;
	Individual* prev_indi2_ptr=NULL;
	pair_id_t prev_pair_id=(pair_id_t)-1;

	assert(size_genotype_row>0);

	LOG( "Read data" );

	/* two pass algorithm: the first algo will store and count the pairs
	 * the second pass will store the data in linkage
	 */
	//for(int side=0; side<2;++side)
		//{
		//LOG("side="<<side);
		for(set<string>::iterator rf=datafiles.begin();
			rf!=datafiles.end();
			++rf)
			{
			//LOG( "Reading data file " << (*rf) << " Step=" << side);
			vector<char*> tokens;
			GZReader in(rf->c_str());
			char* line=NULL;
			size_t nLine=0;


			while((line=in.readLine())!=NULL)
				{
				//LOG(side << " "<< line);
				++nLine;
				if(line[0]==0 || line[0]=='#') continue;
				//if(nLine>1000) break;//TODO
				split(tokens,line,"\t");
				if(tokens.size()!=8) THROW("Illegal number of tokens in "<< (*rf) << " line "<< nLine);

				pair_id_t pair_id=(pair_id)-1;
				Family* fam=NULL;
				Individual* indi1=NULL;
				Individual* indi2=NULL;

				//same couple than previous line
				if(  prev_fam.compare(tokens[0])==0 &&
					 prev_indi1.compare(tokens[1])==0 &&
					 prev_indi2.compare(tokens[2])==0)
					{
					fam=prev_fam_ptr;
					indi1=prev_indi1_ptr;
					indi2=prev_indi2_ptr;
					pair_id=prev_pair_id;
					}
				else
					{
					map<string,Family*>::iterator r=name2family.find(tokens[0]);
					if(r==name2family.end())
						{
						THROW("In "<< (*rf) << " line "<< nLine << " undefined family "<< tokens[0]);
						}
					fam= r->second;
					prev_fam_ptr=fam;
					prev_fam.assign(tokens[0]);

					LOG("fam="<<*fam);
					indi1= fam->findIndividualByName(tokens[1]);
					if(indi1==NULL)
						{
						THROW("In "<< (*rf) << " line "<< nLine << " undefined individual[1] "<< tokens[0]<<":"<< tokens[1]);
						}
					prev_indi1_ptr=indi1;
					prev_indi1.assign(tokens[1]);

					LOG("indi1="<<*indi1);
					indi2= fam->findIndividualByName(tokens[2]);
					if(indi2==NULL)
						{
						THROW("In "<< (*rf) << " line "<< nLine << " undefined individual[2] "<< tokens[0]<<":"<< tokens[2]);
						}
					LOG("indi2="<<*indi2);
					prev_indi2_ptr=indi2;
					prev_indi2.assign(tokens[2]);


					//find the key-ID for this couple
					PairOfIndividual couple(indi1->id,indi2->id);
					Entry keyE((void*)&couple,sizeof(PairOfIndividual));
					auto_ptr<Entry> pairIdEntry=pairOfIndividualsToPairId->get(txn,&keyE,0);
					if(pairIdEntry.get()==NULL)
						{
						THROW("Cannot find pair id for "<< couple);
						}

					//OK, found
					memcpy((void*)&pair_id,pairIdEntry->data(),sizeof(pair_id_t));
					LOG("OK pair id is " << pair_id);
					if(pair_id<0 || pair_id>= this->count_pairs)
						{
						THROW("Runtime error "<<pair_id << " vs " << this->count_pairs);
						}

					prev_pair_id=pair_id;
					}


				/* check chromosome */
				tid_t chromosome_id;
				map<string,tid_t>::iterator rc=chrom2tid.find(tokens[3]);
				if(rc==chrom2tid.end())
					{
					THROW("Runtime error. Undefined chromosome");
					}
				else
					{
					chromosome_id=rc->second;
					LOG("chrom is "<< (int)chromosome_id);
					}



				//if(side==1)
					{
					Entry rsNameE(tokens[4]);
					auto_ptr<Entry> posEntry=this->rsToPositionDb->get(txn,&rsNameE);
					if(posEntry.get()==NULL)
						{
						LOG( "Cannot find position for rs "<< tokens[4] );
						}
					else
						{
						Position position;
						memcpy((void*)&position,posEntry->data(),sizeof(Position));
						Entry positionE((void*)&position,sizeof(Position));
						auto_ptr<Entry> rowEntry= linkageDB->get(txn,&positionE,0);

						if(rowEntry.get()==NULL)
							{
							LOG("new row for "<< tokens[4]);
							for(uint32_t i=0;i< this->count_pairs;++i)
								{
								for(int j=0;j< COUNT_PVALUE;++j)
									{
									genotypeRow[i].pvalue[j]=-1.0f;
									}
								}
							LOG("end new row");
							}
						else
							{
							LOG("row already exist "<< tokens[4]);
							if(rowEntry->size()!=size_genotype_row) THROW("runtime error");
							std::memcpy((void*)genotypeRow,rowEntry->data(),size_genotype_row);

							for(uint32_t i=0;i< this->count_pairs;++i)
									{
									for(int j=0;j< COUNT_PVALUE;++j)
										{
										LOG(i<<"," << j << ":" << genotypeRow[i].pvalue[j]);
										}
									}

							}
						for(int j=0;j< COUNT_PVALUE;++j)
							{
							pvalue_t pvalue=(pvalue_t)atof(tokens[5+j]);
							if( genotypeRow[pair_id].pvalue[j]!=-1.0f &&
								fabs(genotypeRow[pair_id].pvalue[j]-pvalue)> std::numeric_limits<pvalue_t>::epsilon() &&
								max_warning_count>0
								)
								{
								cerr << "Warning duplicate data!! for\n "
									<< 	tokens[4]<<"\n "
									<< position << "\n"
									<< (*indi1) << "\n "
									<< (*indi2)<< "\n "
									<< "pair-id=" << pair_id
									<< " was ["<<j<<"]=" << genotypeRow[pair_id].pvalue[j]
									<< " and now "<< pvalue
									<<endl;
								max_warning_count--;
								if(max_warning_count==0)
									{
									for(uint32_t i=0;i< count_pairs;++i)
										{
										cerr << (i==pair_id?"*":" ")<<"["<< i << "]" << genotypeRow[i] << endl;
										}
									}
								}
							genotypeRow[pair_id].pvalue[j]=pvalue;
							}
						Entry dataL((void*)genotypeRow,size_genotype_row);
						linkageDB->put(txn,&positionE,&dataL,0);
						}
					}


				}
			in.close();

			LOG( (*rf) << ":"<< nLine);
			}
		//LOG("END of side=="<< side );
		//}

	free(genotypeRow);
	genotypeRow=NULL;

	linkageDB->sync();

	}
#endif


void IBDEnv::readData(std::set<std::string>& datafiles)
	{
	int max_warning_count=10;
	/** row of count_pair*Genotype */
	size_t size_genotype_row=0UL;
	Genotype* genotypeRow=NULL;

	string prev_fam;
	Family* prev_fam_ptr=NULL;
	string prev_indi1;
	Individual* prev_indi1_ptr=NULL;
	string prev_indi2;
	Individual* prev_indi2_ptr=NULL;
	pair_id_t prev_pair_id=(pair_id_t)-1;

	LOG( "Read data" );
	/* Initialize the count of pairs */
	this->count_pairs=0;
	/* two pass algorithm: the first algo will store and count the pairs
	 * the second pass will store the data in linkage
	 */
	for(int side=0; side<2;++side)
		{
		LOG("side="<<side);
		if(side==1)
			{
			LOG("alloc row  " << count_pairs);
#ifdef PARANOID
			auto_ptr<Cursor> c= pairOfIndividualsToPairId->cursor();
			while(c->next())
				{
				PairOfIndividual p;
				memcpy(&p,c->key()->data,sizeof(PairOfIndividual));
				LOG(p);
				pair_id_t q;
				memcpy(&q,c->data()->data,sizeof(pair_id_t));
				LOG(q);
				}
			c.reset();


			c= pairIdToPairOfIndividualsDb->cursor();
			while(c->next())
				{
				pair_id_t q;
				memcpy(&q,c->key()->data,sizeof(pair_id_t));
				LOG(q);
				PairOfIndividual p;
				memcpy(&p,c->data()->data,sizeof(PairOfIndividual));
				LOG(p);

				}
			c.reset();
#endif

			/* now we know the number of pairs */
			size_genotype_row=count_pairs*sizeof(Genotype);
			genotypeRow=(Genotype*)safeCalloc(count_pairs,sizeof(Genotype));
			}
		for(set<string>::iterator rf=datafiles.begin();
			rf!=datafiles.end();
			++rf)
			{
			LOG( "Reading data file " << (*rf) << " Step=" << side);
			vector<char*> tokens;
			GZReader in(rf->c_str());
			char* line=NULL;
			size_t nLine=0;

			const int FAMILY1_COLUMN=0;
			const int INDI1_COLUMN=1;
			const int INDI2_COLUMN=3;
			const int CHROM_COLUMN=4;
			const int RS_COLUMN=5;
			const int IBD0_COLUMN=6;
			while((line=in.readLine())!=NULL)
				{
				LOG(side << " "<< line);
				++nLine;
				if(line[0]==0 || line[0]=='#') continue;
				//if(nLine>1000) break;//TODO
				while(line[0]!=0 && line[strlen(line)-1]=='\t')
					{
					line[strlen(line)-1]=0;
					}
				split(tokens,line,"\t");
				if(tokens.size()!=9) THROW("Illegal number of tokens in "<< (*rf) << " line "<< nLine);

				pair_id_t pair_id=(pair_id)-1;
				Family* fam=NULL;
				Individual* indi1=NULL;
				Individual* indi2=NULL;

				//same couple than previous line
				if(  prev_fam.compare(tokens[FAMILY1_COLUMN])==0 &&
					 prev_indi1.compare(tokens[INDI1_COLUMN])==0 &&
					 prev_indi2.compare(tokens[INDI2_COLUMN])==0)
					{
					fam=prev_fam_ptr;
					indi1=prev_indi1_ptr;
					indi2=prev_indi2_ptr;
					pair_id=prev_pair_id;
					}
				else
					{
					map<string,Family*>::iterator r=name2family.find(tokens[FAMILY1_COLUMN]);
					if(r==name2family.end())
						{
						THROW("In "<< (*rf) << " line "<< nLine << " undefined family "<< tokens[FAMILY1_COLUMN]);
						}
					fam= r->second;
					prev_fam_ptr=fam;
					prev_fam.assign(tokens[FAMILY1_COLUMN]);

					LOG("fam="<<*fam);
					indi1= fam->findIndividualByName(tokens[INDI1_COLUMN]);
					if(indi1==NULL)
						{
						THROW("In "<< (*rf) << " line "<< nLine << " undefined individual[1] "<< tokens[FAMILY1_COLUMN]<<":"<< tokens[INDI1_COLUMN]);
						}
					prev_indi1_ptr=indi1;
					prev_indi1.assign(tokens[INDI1_COLUMN]);

					LOG("indi1="<<*indi1);
					indi2= fam->findIndividualByName(tokens[INDI2_COLUMN]);
					if(indi2==NULL)
						{
						THROW("In "<< (*rf) << " line "<< nLine << " undefined individual[2] "<< tokens[FAMILY1_COLUMN]<<":"<< tokens[INDI2_COLUMN]);
						}
					LOG("indi2="<<*indi2);
					prev_indi2_ptr=indi2;
					prev_indi2.assign(tokens[INDI2_COLUMN]);


					//find the key-ID for this couple
					PairOfIndividual couple(indi1->id,indi2->id);
					Entry keyE((void*)&couple,sizeof(PairOfIndividual));
					auto_ptr<Entry> pairIdEntry=pairOfIndividualsToPairId->get(txn,&keyE,0);
					if(pairIdEntry.get()!=NULL)
						{
						//OK, found
						memcpy((void*)&pair_id,pairIdEntry->data(),sizeof(pair_id_t));
						LOG("OK pair id is " << pair_id);
						if(side==1 && pair_id>= this->count_pairs)
							{
							THROW("Runtime error "<<pair_id << " vs " << this->count_pairs);
							}


						}
					else
						{
						LOG("Cannot find pair id for "<< couple);
						if(side==1) THROW("Runtime error cannot find pair id for  "<< couple);
						pair_id=this->count_pairs;
						this->count_pairs++;
						if(this->count_pairs%100==0) LOG("pairs inserted "<< count_pairs);
						//update databases
						LOG("inserting new pair "<< pair_id);
						Entry dataE((void*)&pair_id,sizeof(pair_id_t));
						pairIdToPairOfIndividualsDb->put(txn,&dataE,&keyE,0);
						pairOfIndividualsToPairId->put(txn,&keyE,&dataE,0);
						assert(pairIdToPairOfIndividualsDb->contains(txn,&dataE));
						assert(pairOfIndividualsToPairId->contains(txn,&keyE));
						}
					prev_pair_id=pair_id;
					}


				/* check chromosome */
				tid_t chromosome_id;
				map<string,tid_t>::iterator rc=chrom2tid.find(tokens[CHROM_COLUMN]);
				if(rc==chrom2tid.end())
					{
					THROW("Runtime error. Undefined chromosome");
					}
				else
					{
					chromosome_id=rc->second;
					LOG("chrom is "<< (int)chromosome_id);
					}



				if(side==1)
					{
					Entry rsNameE(tokens[RS_COLUMN]);
					auto_ptr<Entry> posEntry=this->rsToPositionDb->get(txn,&rsNameE);
					if(posEntry.get()==NULL)
						{
						LOG( "Cannot find position for rs "<< tokens[RS_COLUMN] );
						}
					else
						{
						Position position;
						memcpy((void*)&position,posEntry->data(),sizeof(Position));
						Entry positionE((void*)&position,sizeof(Position));
						auto_ptr<Entry> rowEntry= linkageDB->get(txn,&positionE,0);

						if(rowEntry.get()==NULL)
							{
							LOG("new row for "<< tokens[RS_COLUMN]);
							for(uint32_t i=0;i< this->count_pairs;++i)
								{
								for(int j=0;j< COUNT_PVALUE;++j)
									{
									genotypeRow[i].pvalue[j]=-1.0f;
									}
								}
							LOG("end new row");
							}
						else
							{
							LOG("row already exist "<< tokens[RS_COLUMN]);
							if(rowEntry->size()!=size_genotype_row) THROW("runtime error");
							std::memcpy((void*)genotypeRow,rowEntry->data(),size_genotype_row);

							for(uint32_t i=0;i< this->count_pairs;++i)
									{
									for(int j=0;j< COUNT_PVALUE;++j)
										{
										LOG(i<<"," << j << ":" << genotypeRow[i].pvalue[j]);
										}
									}

							}
						for(int j=0;j< COUNT_PVALUE;++j)
							{
							pvalue_t pvalue=(pvalue_t)atof(tokens[IBD0_COLUMN+j]);
							if( genotypeRow[pair_id].pvalue[j]!=-1.0f &&
								fabs(genotypeRow[pair_id].pvalue[j]-pvalue)> std::numeric_limits<pvalue_t>::epsilon() &&
								max_warning_count>0
								)
								{
								cerr << "Warning duplicate data!! for\n "
									<< 	tokens[RS_COLUMN]<<"\n "
									<< position << "\n"
									<< (*indi1) << "\n "
									<< (*indi2)<< "\n "
									<< "pair-id=" << pair_id
									<< " was ["<<j<<"]=" << genotypeRow[pair_id].pvalue[j]
									<< " and now "<< pvalue
									<<endl;
								max_warning_count--;
								if(max_warning_count==0)
									{
									for(uint32_t i=0;i< count_pairs;++i)
										{
										cerr << (i==pair_id?"*":" ")<<"["<< i << "]" << genotypeRow[i] << endl;
										}
									}
								}
							genotypeRow[pair_id].pvalue[j]=pvalue;
							}
						Entry dataL((void*)genotypeRow,size_genotype_row);
						linkageDB->put(txn,&positionE,&dataL,0);
						}
					}


				}
			in.close();

			LOG( (*rf) << ":"<< nLine);
			}
		if(side==1)
			{
			free(genotypeRow);
			genotypeRow=NULL;
			}
		LOG("END of side=="<< side );
		}

	linkageDB->sync();
	pairIdToPairOfIndividualsDb->sync();
	pairOfIndividualsToPairId->sync();
	}



Application::Application():db_home(getenv("DB_HOME"))
	{
	}

Application::~Application()
	{
	ibdEnv.close();
	}

IBDEnv* Application::env()
	{
	return &ibdEnv;
	}

int Application::checkHome()
	{
	if(db_home==NULL)
			{
			cerr << "--db-home undefined" << endl;
			return EXIT_FAILURE;
			}
	struct stat sb;
	if(stat(db_home, &sb) != 0)
		{
		cerr << db_home << " doesn't exist.\n";
		return EXIT_FAILURE;
		}
	if(!S_ISDIR(sb.st_mode))
		{
		cerr << db_home << " is not a directory.\n";
		return EXIT_FAILURE;
		}
	return EXIT_SUCCESS;
	}

ApplicationBuilder::ApplicationBuilder():
		reskinfile(NULL),
		pedfile(NULL),
		mapfile(NULL),
		datafiles(NULL),
		bedfile(NULL)
	{

	}

ApplicationBuilder::~ApplicationBuilder()
	{

	}

void ApplicationBuilder::usage(std::ostream& out)
	{
	out << "options\n" <<
		" --db-home (berkeydb home) default $DB_HOME " << (db_home==NULL?"":db_home)<< " \n" <<
		" --ped <pedigree file> \n"<<
		" --reskin <reskin file> \n"<<
		" --map <map file (Axiom_GW_Hu_SNP)>\n"<<
		"  OR\n"<<
		" --bed <BED>  \n"<<
		" --data <file-containing-the-path-to-the-data-files> \n"<<
		std::endl
		;
	}

int ApplicationBuilder::main(int argc,char** argv)
	{
	if(argc==1)
		{
		usage(cout);
		return EXIT_SUCCESS;
		}
	int optind=1;
	while(optind < argc)
		{

		if(std::strcmp(argv[optind],"-h")==0)
			{
			usage(cout);
			std::exit(EXIT_FAILURE);
			}
		else if(std::strcmp(argv[optind],"--db-home")==0 && optind+1<argc)
			{
			db_home=argv[++optind];
			}
		else if(std::strcmp(argv[optind],"--ped")==0 && optind+1<argc)
			{
			pedfile=argv[++optind];
			}
		else if(std::strcmp(argv[optind],"--data")==0 && optind+1<argc)
			{
			datafiles=argv[++optind];
			}
		else if(std::strcmp(argv[optind],"--map")==0 && optind+1<argc)
			{
			mapfile=argv[++optind];
			}
		else if(std::strcmp(argv[optind],"--bed")==0 && optind+1<argc)
			{
			bedfile=argv[++optind];
			}
		else if(std::strcmp(argv[optind],"--reskin")==0 && optind+1<argc)
			{
			reskinfile=argv[++optind];
			}
		else if(std::strcmp(argv[optind],"--debug")==0)
			{
			ibdEnv.debug=true;
			}
		else if(std::strcmp(argv[optind],"--")==0)
			{
			++optind;
			break;
			}
		else if(argv[optind][0]=='-')
			{
			std::cerr << "unknown option " << argv[optind] << std::endl;
			std::exit(EXIT_FAILURE);
			}
		else
			{
			break;
			}
		++optind;
		}
	if(checkHome()!=0) return EXIT_FAILURE;
	if(pedfile==NULL)
		{
		cerr << "undefined pedigree file." << endl;
		return EXIT_FAILURE;
		}
	if(mapfile==NULL && bedfile==NULL)
		{
		cerr << "undefined map or bed file." << endl;
		return EXIT_FAILURE;
		}
	if(reskinfile==NULL)
		{
		cerr << "undefined reskin file." << endl;
		return EXIT_FAILURE;
		}
	set<string> datafilesset;
	if(datafiles!=NULL)
		{
		GZReader r(datafiles);
		char* line;
		while((line=r.readLine())!=NULL)
			{
			if(line[0]==0 || line[0]=='#') continue;
			datafilesset.insert(line);
			}
		r.close();
		}
	if(datafilesset.empty())
		{
		cerr << "undefined data file." << endl;
		return EXIT_FAILURE;
		}

	ibdEnv.open(db_home,false);
	try
		{
		ibdEnv.readPedigree(pedfile);
		ibdEnv.readData(datafilesset);
		if(mapfile!=0)
		    {
		    ibdEnv.readMapping(mapfile);
		    }
		else
		    {
		    ibdEnv.readBed(bedfile);
		    }
		ibdEnv.readReskin(reskinfile);
		}
	catch(std::exception& err)
		{
		cerr << "An error occured:" << err.what()<<  endl;
		}
	catch(...)
		{
		cerr << "An error occured" << endl;
		}
	ibdEnv.close();
	return EXIT_SUCCESS;
	}



ApplicationReader::ApplicationReader():
		output_mode(ONE_PAIR_PER_LINE),
		print_self_self(true),
		print_theoritical_proba(true),
		min_pvalue_ibd0(-1.0),
		output(&cout),
		errstream(&cerr)
	{
	for(int i=0;i<COUNT_PVALUE;++i )
	    {
	    accept_pvalue_at[i]=true;
	    }
	}

ApplicationReader::~ApplicationReader()
	{

	}


#define TRACE cout << __FILE__ << ":"<< __LINE__ << ":" << endl;cerr << __FILE__ << ":"<< __LINE__ << ":" << endl

#define PRINTERROR(a) if(this->errstream!=NULL) do { (*errstream) << __FILE__ << ":" << __LINE__ << ":" << a << std::endl;} while(0)
void ApplicationReader::scanRegion(GenomicPosition* genomic)
    {

    if(genomic->chromStart>genomic->chromEnd) return;

    map<string,tid_t>::iterator r=ibdEnv.chrom2tid.find(genomic->chromosome);
    if(r==ibdEnv.chrom2tid.end())
		{

    	PRINTERROR("unknown chromosome \""<< genomic->chromosome << "\"" );
		return;
		}
    Position position;
    position.chrom=r->second;
    position.pos=genomic->chromStart;
    Entry key((void*)&position,sizeof(Position));
    auto_ptr<Cursor> c=ibdEnv.pos2markerDb->cursor(NULL,(u_int32_t)0);
    bool first=true;
    for(;;)
		{

		if(first)
			{
			first=false;
			DBT dbtK;
			DBT dbtV;
			memset(&dbtK,0,sizeof(DBT));
			memset(&dbtV,0,sizeof(DBT));
			dbtK.data=(void*)&position;
			dbtK.size=sizeof(Position);
			if(c->get(&dbtK,&dbtV,DB_SET_RANGE)!=0)
				{
				break;
				}
			}
		else
			{
			if(!c->next())
				{
				break;
				}
			}

		Position curr;
		memcpy((void*)&curr,c->key()->data,sizeof(Position));



		if(curr.chrom!=position.chrom ||
		   curr.pos> genomic->chromEnd)
			{

			break;
			}

		scanPosition(&curr);

		}
    }

void ApplicationReader::resultHeader1()
    {
    std::ostream& out=(*output);
    switch(output_mode)
		{
		case ALL_PAIRS_ON_LINE:
			{
			out << "#rs\t";
			out << "chrom"<< "\t";
			out << "pos" << "\t";
			out << "geneticMaps" ;
			auto_ptr<Cursor> c=ibdEnv.pairIdToPairOfIndividualsDb->cursor(ibdEnv.txn,(uint32_t)0);
			while(c->next())
				{
				pair_id_t pair_id;
				c->keyCopyTo((void*)&pair_id,0,sizeof(pair_id_t));
				if(!accept_pair_id(pair_id)) continue;

				PairOfIndividual couple;
				c->dataCopyTo((void*)&couple,0,sizeof(PairOfIndividual));

				auto_ptr<Genotype> theoritical= ibdEnv.findTheoriticalIBDByPairId(pair_id);

				CoupleOfIndividual p=ibdEnv.findIndividualsByPair(&couple);
				if(p.first==NULL || p.second==NULL) THROW("Bad state "<< pair_id);
				if(!accept_individual_id(p.first->id)) continue;
				if(!accept_individual_id(p.second->id)) continue;
				if(!accept_family_id(p.first->family->id)) continue;
				if(!accept_pair_self_self() && couple.first==couple.second) continue;

				for(int i=0;i< COUNT_PVALUE;++i)
					{
					if(!accept_pvalue_at_index(i)) continue;

						out <<"\t";
						out << p.first->family->name;
						out << "|";
						out << p.first->name;
						out << "|";
						out << p.second->name;
						if(print_theoritical_proba && theoritical.get()!=NULL)
							{
							out << "|theoIBD"<<i<<":"<< theoritical->pvalue[i];
							}


					}
				}
			c.reset();
			if(min_pvalue_ibd0>0)
				{
				out << "\tcount(IBD0>"<< min_pvalue_ibd0<<")\tcount(IBD0>"<< min_pvalue_ibd0<<")/pairs";
				}
			out << endl;
			break;
			}
		case ONE_PAIR_PER_LINE:
			out << "#rs\t";
			out << "chrom"<< "\t";
			out << "pos" << "\t";
			out << "geneticMaps" << "\t";
			out << "family" << "\t";
			out << "indi1.name" << "\t";
			out << "indi1.father" << "\t";
			out << "indi1.mother" << "\t";
			out << "indi1.gender" << "\t";
			out << "indi1.status" << "\t";
			out << "indi2.name" << "\t";
			out << "indi2.father" << "\t";
			out << "indi2.mother" << "\t";
			out << "indi2.gender" << "\t";
			out << "indi2.status";
			for(int i=0;i< COUNT_PVALUE;++i)
				{
				if(!accept_pvalue_at_index(i)) continue;
				if(print_theoritical_pvalue())
						{
						out << "\ttheoIBD"<< i;
						}
				out << "\texpIBD"<< i;
				}
			out << endl;
			break;
		default: THROW("Illegal state "<< output_mode);break;
		}
    }

void ApplicationReader::readRsList(Reader* reader)
	{
	char* line=NULL;

	while((line=reader->readLine())!=NULL)
	    {
	    if(line[0]==0 || line[0]=='#') continue;

	    Entry rsIDE(line);

	    /* position */
	    auto_ptr<Entry> posEntry= ibdEnv.rsToPositionDb->get(ibdEnv.txn,&rsIDE,0);
	    if(posEntry.get()==NULL)
			{
	    		PRINTERROR( "Cannot find " << line << " in database.");
			continue;
			}

	    Position position;
	    memcpy((void*)&position,posEntry->data(),sizeof(Position));

	    scanPosition(&position);
	    }
	}

bool ApplicationReader::print_theoritical_pvalue()
		{
		return this->print_theoritical_proba;
		}

bool ApplicationReader::accept_family_id(family_id_t famid)
	{
	return !restric_family_ids.empty() &&
	    restric_family_ids.find(famid)==restric_family_ids.end()
    		? false: true;
	}

bool ApplicationReader::accept_individual_id(individual_id_t indid)
	{
	return !restric_individuals_ids.empty() &&
		restric_individuals_ids.find(indid)==restric_individuals_ids.end()
    		? false: true;
	}

bool ApplicationReader::accept_pair_id(pair_id_t pairid)
	{
	return !restric_pair_ids.empty() &&
		restric_pair_ids.find(pairid)==restric_pair_ids.end()
		? false: true;
	}

bool ApplicationReader::accept_marker(const Marker* marker)
	{
	return true;
	}

bool ApplicationReader::accept_position(const Position* pos)
	{
	return true;
	}

bool ApplicationReader::accept_pvalue_at_index(int p)
    {
    assert(p>=0 && p< COUNT_PVALUE);
    return accept_pvalue_at[p];
    }

bool ApplicationReader::accept_pair_self_self()
    {
    return print_self_self;
    }

#define ECHO_MARKER_INFO out << marker->rsId<< "\t";\
	out << rchrom->second << "\t";\
	out << position->pos << "\t";\
	for(uint32_t i=0;i< marker->geneticMaps.size();++i)\
		{\
		if(i>0) out << "|";\
		out <<  marker->geneticMaps.at(i);\
		}


int ApplicationReader::scanPosition(const Position* position)
    {
    std::ostream& out=(*output);

    if(!accept_position(position))
	    {
	    return EXIT_SUCCESS;
	    }

    /** position */
    map<tid_t,string>::iterator rchrom=ibdEnv.id2chromosome.find(position->chrom);
    if(rchrom==ibdEnv.id2chromosome.end())
		{
    		PRINTERROR("Bad tid " << position->chrom << ".");
		return EXIT_FAILURE;
		}

    Entry posEntry((void*)position,sizeof(Position));
    /* marker */
    auto_ptr<Entry> markerEntry= ibdEnv.pos2markerDb->get(ibdEnv.txn,&posEntry,0);
    if(markerEntry.get()==NULL)
		{
    		PRINTERROR("Cannot find description for marker at " << *position);
		return EXIT_FAILURE;
		}


    auto_ptr<Marker> marker(ibdEnv.entryToMarker(markerEntry.get()));

    if(!accept_marker(marker.get()))
        	{
        	return EXIT_SUCCESS;
        	}

    if(memcmp((void*)&(marker->position),(void*)position,sizeof(Position))!=0)
		{
		THROW("Illegal state "<< *position);
		}


    /* linkage */

    auto_ptr<Entry> linkageRowEntry= ibdEnv.linkageDB->get(ibdEnv.txn,&(posEntry),0);

    if(linkageRowEntry.get()==NULL)
		{
    		PRINTERROR("position "<< (*position) << " is not in IBD data.");
		return EXIT_FAILURE;
		}



    if(linkageRowEntry->size()< ibdEnv.count_pairs*sizeof(Genotype))
		{
		THROW("Illegal state "<< linkageRowEntry->size());
		}

    switch(output_mode)
		{
		case ALL_PAIRS_ON_LINE:
			{
			ECHO_MARKER_INFO
			break;
			}
		default:break;
		}
    int count_having_ibd0_gt=0;
    int count_pairs_printed=0;
    pair_id_t pair_id=0;
    auto_ptr<Cursor> c=ibdEnv.pairIdToPairOfIndividualsDb->cursor(ibdEnv.txn,(uint32_t)0);
    const Genotype* genotypeRow=(const Genotype*)linkageRowEntry->data();
    while(c->next())
		{

    	//check pair id are sorted
		assert(sizeof(pair_id_t)==c->key()->size);
		if(memcmp((void*)&pair_id,c->key()->data,c->key()->size)!=0)
			{
			pair_id_t p2;
			memcpy((void*)&p2,c->key()->data,c->key()->size);

			THROW("Illegal state expected pair_id="<< pair_id << " vs found cursor:" <<  p2 );
			//++pair_id;
			}

		if(!accept_pair_id(pair_id))
			{
			++pair_id;
			continue;
			}



		PairOfIndividual couple;
		assert(sizeof(PairOfIndividual)==c->data()->size);
		memcpy((void*)&couple,c->data()->data,c->data()->size);

		if(!accept_pair_self_self() &&
		    couple.first == couple.second  )
		    {
		    ++pair_id;
		    continue;
		    }

		if(!accept_individual_id(couple.first))
			{
			++pair_id;
			continue;
			}

		if(!accept_individual_id(couple.second))
			{
			++pair_id;
			continue;
			}

		IndividualPtr indi1=ibdEnv.findIndividualById(couple.first);
		if(indi1==NULL)
			{
			THROW("Cannot find individual "<< couple.first);
			}



		//do we need to scan his family ?
		if(!accept_family_id(indi1->family->id))
			{
			++pair_id;
			continue;
			}

		Individual* indi2=ibdEnv.findIndividualById(couple.second);
		if(indi2==NULL)
		    {
		    THROW("Cannot find individual "<< couple.second);
		    }

		Genotype genotype;

		memcpy((void*)&genotype,(void*)(&genotypeRow[pair_id]),sizeof(Genotype));




		switch(output_mode)
		    {
		    case ALL_PAIRS_ON_LINE:
				{
				for(int i=0;i< COUNT_PVALUE;++i)
					{
					if(!accept_pvalue_at_index(i)) continue;
					out << "\t";
					out << genotype.pvalue[i];
					}
				if(min_pvalue_ibd0>0)
					{
					count_pairs_printed++;
					if(genotype.ibd0()> min_pvalue_ibd0)
						{
						count_having_ibd0_gt++;
						}
					}
				break;
				}
		    case ONE_PAIR_PER_LINE:
				{
				bool hasData=false;
				for(int n_p=0;n_p< COUNT_PVALUE;++n_p)
					{
					if(!accept_pvalue_at_index(n_p)) continue;
					if(genotype.pvalue[n_p]>=0)
						{
						hasData=true;
						break;
						}
					}

				if(hasData)
					{
					auto_ptr<Genotype> theoritical= ibdEnv.findTheoriticalIBDByPairId(pair_id);
					ECHO_MARKER_INFO
					out << "\t";
					out << indi1->family->name << "\t";
					out << indi1->name << "\t";
					out << indi1->father << "\t";
					out << indi1->mother << "\t";
					out << indi1->gender << "\t";
					out << indi1->status << "\t";
					out << indi2->name << "\t";
					out << indi2->father << "\t";
					out << indi2->mother << "\t";
					out << indi2->gender << "\t";
					out << indi2->status ;
					for(int i=0;i< COUNT_PVALUE;++i)
						{
						if(!accept_pvalue_at_index(i)) continue;
						if(print_theoritical_pvalue())
							{
							out << "\t";
							if(theoritical.get()==NULL)
								{
								out << "N/A";
								}
							else
								{
								out << theoritical->pvalue[i];
								}
							}
						out << "\t";
						out << genotype.pvalue[i];
						}
					out << endl;
					}
				break;
				}
		    default:
		    	{
		    	THROW("Uhh??");
		    	break;
		    	}
		    }


		pair_id++;
		}//end of loop over pair
    c->close();

    switch(output_mode)
    	{
    	case ALL_PAIRS_ON_LINE:
    	    {
    	    if(min_pvalue_ibd0>0)
						{
    	    			out << "\t";
						out << count_having_ibd0_gt;
						out << "\t";
						out << ((double)count_having_ibd0_gt/(double)count_pairs_printed);
						}
    	    out << endl;
    	    break;
    	    }
    	default:break;
    	}
    return EXIT_SUCCESS;
    }

void ApplicationReader::proba_christian()
	{


	/** loop over all the markers */
	auto_ptr<Cursor> cpos=ibdEnv.linkageDB->cursor(ibdEnv.txn,(uint32_t)0);
	while(cpos->next())
		{
		vector<CoupleOfIndividual> pairsAffected;
		vector<CoupleOfIndividual> pairsUnaffected;
		map<pvalue_t,uint32_t> nbTypePairesAA;
		map<pvalue_t,uint32_t> nbTypePairesAU;
		/** loop over all the pairs */
		auto_ptr<Cursor> c=ibdEnv.pairIdToPairOfIndividualsDb->cursor(ibdEnv.txn,(uint32_t)0);
		while(c->next())
				{
				pair_id_t pair_id=0;
				memcmp((void*)&pair_id,c->key()->data,c->key()->size);
				PairOfIndividual couple;
				memcpy((void*)&couple,c->data()->data,c->data()->size);
				CoupleOfIndividual pairIndis=ibdEnv.findIndividualsByPair(&couple);
				auto_ptr<Genotype> genotype=ibdEnv.findTheoriticalIBDByPairId(pair_id);

				switch(pairIndis.concordant())
					{
					case 2:
						{
						++nbTypePairesAA[genotype->ibd0()];
						pairsAffected.push_back(pairIndis);
						break;
						}
					case 1:
						{
						++nbTypePairesAU[genotype->ibd0()];
						pairsUnaffected.push_back(pairIndis);
						break;
						}
					default:break;
					}

				}
		c.reset();

		}
	cpos.reset();
	}


void ApplicationReader::usage(std::ostream& out)
	{
	out << "options\n" <<
			" --db-home (berkeydb home) default $DB_HOME " << (db_home==NULL?"":db_home)<< " \n" <<
			" --position chrom:start-end\n"
			" --all dump all\n"
			" -o <output-file> (or stdout if omitted).\n"
			" --restrict-family <family-name> (optional, can be called several times).\n"
			" --restrict-individual  <family-name> <individual-name> (optional, can be called several times).\n"
			" --no-self don't print self-self pairs.\n"
			" --no-theo don't print theoritical proba.\n"
			" --mode '0' : one pair per line (=default) '1' all pairs on one line.\n"
			" --minibd0 (value) print a count of the pair having IBD0>value.\n";
	for(int i=0;i< COUNT_PVALUE;++i)
	    {
	    out << " --pvalue"<<(i+1)<< " disable pvalue["<<(i+1)<< "]" << endl;
	    }
	out <<  "(files with rs name | stdin )\n"<<
			std::endl
			;
	}

int ApplicationReader::main(int argc,char** argv)
	{
	if(argc==1)
		{
		usage(cout);
		return EXIT_SUCCESS;
		}
	bool dump_all=false;
	int optind=1;
	char* position=NULL;
	char* fileout=NULL;
	auto_ptr<fstream> oustream;
	set<string> restrict_family_names;
	vector<string> restrict_individual_names;
	while(optind < argc)
		{

		if(std::strcmp(argv[optind],"-h")==0)
		    {
		    usage(cout);
		    std::exit(EXIT_FAILURE);
		    }
		else if(std::strcmp(argv[optind],"--mode")==0 && optind+1<argc)
		    {
		    switch(atoi(argv[++optind]))
				{
				case 1: this->output_mode= ALL_PAIRS_ON_LINE;break;
				default: this->output_mode= ONE_PAIR_PER_LINE;break;
				}
		    }
		else if(std::strcmp(argv[optind],"--no-self")==0)
		    {
		    print_self_self=false;
		    }
		else if(std::strcmp(argv[optind],"--no-theo")==0)
			{
			print_theoritical_proba=false;
			}
		else if(std::strncmp(argv[optind],"--pvalue",8)==0)
		    {
		    int pv=atoi(&argv[optind][8])-1;
		    if(pv>=0 && pv< COUNT_PVALUE)
				{
				accept_pvalue_at[pv]=false;
				}
		    else
				{
				std::cerr << "unknown option " << argv[optind] << std::endl;
				return EXIT_FAILURE;
				}
		    }
		else if(std::strcmp(argv[optind],"--minibd0")==0 && optind+1<argc)
			{
			min_pvalue_ibd0=atof(argv[++optind]);
			}
		else if(std::strcmp(argv[optind],"--db-home")==0 && optind+1<argc)
		    {
		    db_home=argv[++optind];
		    }
		else if(std::strcmp(argv[optind],"--position")==0 && optind+1<argc)
		    {
		    position=argv[++optind];
		    }
		else if(std::strcmp(argv[optind],"--all")==0 )
			{
			dump_all=true;
			}
		else if(std::strcmp(argv[optind],"-o")==0 && optind+1<argc)
			{
			fileout=argv[++optind];
			}
		else if(std::strcmp(argv[optind],"--restrict-family")==0 && optind+1<argc)
			{
			restrict_family_names.insert(argv[++optind]);
			}
		else if(std::strcmp(argv[optind],"--restrict-individual")==0 && optind+2<argc)//2!
			{
			restrict_individual_names.push_back(argv[++optind]);//fam
			restrict_individual_names.push_back(argv[++optind]);//indi
			}
		else if(std::strcmp(argv[optind],"--debug")==0)
			{
			ibdEnv.debug=true;
			}
		else if(std::strcmp(argv[optind],"--")==0)
			{
			++optind;
			break;
			}
		else if(argv[optind][0]=='-')
			{
			std::cerr << "unknown option " << argv[optind] << std::endl;
			return (EXIT_FAILURE);
			}
		else
			{
			break;
			}
		++optind;
		}
	if(checkHome()!=0) return EXIT_FAILURE;
	if((position!=NULL || dump_all) && optind!=argc)
	    {
	    cerr << "Illegal number of arguments (--position was used)"<< endl;
	    return EXIT_FAILURE;
	    }
	if(fileout!=NULL)
		{
		oustream.reset(new fstream(fileout,ios::out));
		if(!oustream->is_open())
			{
			cerr << "Cannot open \""<< fileout << "\" " << strerror(errno) << endl;
			return EXIT_FAILURE;
			}
		this->output=oustream.get();
		}
	ibdEnv.open(db_home,true);

	for(set<string>::iterator r=restrict_family_names.begin();
			r!=restrict_family_names.end();++r)
		{
		FamilyPtr r2=ibdEnv.findFamilyByName(r->c_str());
		if(r2==NULL)
			{
			cerr << "Unknown family \""<< (*r) << "\"" << endl;
			return EXIT_FAILURE;
			}
		this->restric_family_ids.insert(r2->id);
		}

	for(size_t i=0;i+1< restrict_individual_names.size();i+=2)
		{
		FamilyPtr f2=ibdEnv.findFamilyByName(restrict_individual_names[i].c_str());
		if(f2==NULL)
			{
			cerr << "Unknown family \""<< restrict_individual_names[i] << "\"" << endl;
			return EXIT_FAILURE;
			}
		IndividualPtr indiptr= f2->findIndividualByName(restrict_individual_names[i+1].c_str());
		if(indiptr==NULL)
			{
			cerr << "Unknown individual \""<< restrict_individual_names[i+1]
			     << "\" for family :"
			     << f2->name
			     << endl;
			cerr << "Available:\n";
			for(map<int,IndividualPtr>::iterator r=f2->id2individual.begin();
				r!=f2->id2individual.end();
				++r)
				{
				cerr << " \""<< r->second->name << "\""<<endl;
				}
			for(map<string,IndividualPtr>::iterator r=f2->name2individual.begin();
							r!=f2->name2individual.end();
							++r)
				{
				cerr << "-\""<< r->second->name << "\""<<endl;
				}
			return EXIT_FAILURE;
			}
		this->restric_individuals_ids.insert(indiptr->id);
		}

	try
		{
		if(dump_all)
			{
			resultHeader1();
			for(map<tid_t,string>::iterator r= ibdEnv.id2chromosome.begin();
					r!=ibdEnv.id2chromosome.end();++r)
				{
				GenomicPosition gp((*r).second,0,
						gp.chromEnd=numeric_limits<GenomicPosition::position_t>::max()-100 //avoid side effect
					);
				scanRegion(&gp);
				}
			}
		else if(position!=NULL)
		    {
		    auto_ptr<GenomicPosition> ptr=GenomicPosition::parse(position);
		    if(ptr.get()==NULL)
			    {
			    cerr << "Cannot parse position "<< position << endl;
			    return EXIT_FAILURE;
			    }
		    resultHeader1();
		    scanRegion(ptr.get());
		    }
		else if(optind==argc)
			{
			StdioReader in(stdin);
			resultHeader1();
			readRsList(&in);
			}
		else
			{
			while(optind<argc)
				{
				char* filename=argv[optind++];
				GZReader in(filename);
				resultHeader1();
				readRsList(&in);
				in.close();
				}
			}
		}
	catch(std::exception& err)
		{
		cerr << "An error occured:" << err.what()<<  endl;
		}
	catch(...)
		{
		cerr << "An error occured" << endl;
		}
	ibdEnv.close();
	if(oustream.get()!=NULL)
		{
		oustream->flush();
		oustream->close();
		this->output=&cout;
		oustream.reset();
		}
	return EXIT_SUCCESS;
	}


